#include "yocto/utest/run.hpp"
#include "kernel/solution.hpp"

#include "yocto/ios/ocstream.hpp"
#include "yocto/lua/lua-state.hpp"
#include "yocto/lua/lua-config.hpp"
#include "yocto/code/rand.hpp"

using namespace yocto;


YOCTO_UNIT_TEST_IMPL(dosage)
{
	_rand.wseed();
	
	if( argc <= 1 )
		throw exception("%s: need script.lua", argv[0] );
	Library     lib;
	ChemSys     chs(lib);
	Initializer ini(lib);
	
	Lua::State VM;
	lua_State *L = VM();
	Lua::Config::DoFile( L, argv[1] );
	lib.load(L,"library");
	lib.close();
	
	chs.load(L,"chemsys");
	chs.compile();
	
	ini.load(L,"ini");
	ini.compile();
	
	
	Real     Volume = Lua::Config::Get<lua_Number>( L, "Volume" );
	Real     Cb     = Lua::Config::Get<lua_Number>( L, "Cb"     );
	Real     Vb     = Lua::Config::Get<lua_Number>( L, "Vb"     );
	
	Solution S(lib);
	ini.setup(chs, 0);
	S.load_from( chs );
	S.display(std::cerr);
	std::cerr << "pH=" << S.pH() << std::endl;
	const size_t Np  = 500;
	const Real   dV  = Vb  / Np;
	const Real   dNb = Cb * dV;
	
	ios::ocstream fp( "dosage.dat", false );
	fp("0 %g\n", S.pH() );
	
	for( size_t i=1; i <= Np; ++i )
	{
		const Real Vadd = i*dV;
		const Real Vold = Volume + (i-1) * dV; 
		const Real Vtot = Volume + Vadd;
		const Real dC   = dNb / Vtot;
		//-- dilution
		for( size_t j=1; j <= S.size; ++j )
		{
			Real &Cj = S[j];
			const Real Nj = Cj * Vold;
			Cj = Nj/Vtot;
		}
		
		//-- base
		S[ "Na+" ] += dC;
		S[ "HO-" ] += dC;
		
		//-- reaction
		S.save_into(chs);
		chs.solveGamma(0);
		S.load_from(chs);
		std::cerr << "pH=" << S.pH() << std::endl;
		fp("%g %g\n", Vadd, S.pH() );
	}
	
	
}
YOCTO_UNIT_TEST_DONE()